Preparation

Packages

library(tidyverse) 
library(openxlsx) #read excel file
library(ggrepel) #label text
library(vegan) #Biodiversity Index
library(sf) #Map
library(scales)
library(RColorBrewer)
library(treemapify)

Call custom function

source("Custom_function.R")

Read data

Read Bidoup - Nui Ba data with read.xlsx() function

#library(openxlsx)
filename <- "Plot_NC-BD_20230925_v4.xlsx"
# Bidoup - Nui Ba plot
BDplot <- read.xlsx(filename, sheet = "BD_census2")

We practice with 2ha data, in which BD-1 (1ha) is concentrated (100x100m) and BD-2 (1ha) is taken from 25 scattered subplots.

Community (Species composition, size-class structure and diversity)

Number of family, genus and species

Use n_distinct() counts the number of distinct value in column: family, genus, scientificName. And group_by() to count at each plot.

BDplot |>  
  group_by(PlotID) |> 
  summarise("family" = n_distinct(family),
            "genus" = n_distinct(genus),
            "species" = n_distinct(scientificName),
            "stem" = n())
## # A tibble: 2 × 5
##   PlotID family genus species  stem
##   <chr>   <int> <int>   <int> <int>
## 1 BD-1       46    70     124 12207
## 2 BD-2       53    94     169 11575

Forest structure

We usually group trees into 3 size-class of DBH:

  • small tree [1-5)cm;

  • medium tree [5-20)cm;

  • big tree >=20 cm

BDplot <- BDplot |>
  mutate(
    dbh_cm = dbh / 10,
    # dbh value in mm
    gDBH = case_when(dbh_cm >= 20 ~ ">=20", 
                     dbh_cm >= 5 ~ "[5-20)", 
                     TRUE ~ "[1-5)")
  )
head(BDplot)

We count the number of stems at each diameter class and each plot

(stemByDBH <- BDplot |>
   group_by(PlotID, gDBH) |> 
   summarise(n = n()) |> 
    mutate(freq = round(n / sum(n), 3)) |>
    arrange(desc(n)))
## `summarise()` has grouped output by 'PlotID'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 4
## # Groups:   PlotID [2]
##   PlotID gDBH       n  freq
##   <chr>  <chr>  <int> <dbl>
## 1 BD-1   [1-5)   9736 0.798
## 2 BD-2   [1-5)   9478 0.819
## 3 BD-1   [5-20)  2057 0.169
## 4 BD-2   [5-20)  1749 0.151
## 5 BD-1   >=20     414 0.034
## 6 BD-2   >=20     348 0.03

Figure

ggplot(BDplot, aes(factor(gDBH), fill = factor(PlotID))) +
  geom_bar(position = "dodge") +
  theme_bw(14) +
  labs(x = "DBH (cm)", y = "Nb of stem", fill = "PlotID")

We can also count the number of species at each diameter class.

(
  SpByDBH <- BDplot  |>
    group_by(PlotID, gDBH) |>
    summarise(
      "family" = n_distinct(family),
      "genus" = n_distinct(genus),
      "species" = n_distinct(scientificName)
    )
)
## `summarise()` has grouped output by 'PlotID'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 5
## # Groups:   PlotID [2]
##   PlotID gDBH   family genus species
##   <chr>  <chr>   <int> <int>   <int>
## 1 BD-1   >=20       28    35      53
## 2 BD-1   [1-5)      45    69     120
## 3 BD-1   [5-20)     42    61      96
## 4 BD-2   >=20       33    45      64
## 5 BD-2   [1-5)      53    91     160
## 6 BD-2   [5-20)     48    74     121

Figure

# Convert to longer table
SpByDBH <- SpByDBH |>
  tidyr::pivot_longer(-c(PlotID,gDBH), names_to = "rank", values_to = "n")

Plot BD-1

SpByDBH |> filter(PlotID=="BD-1") |> 
  ggplot(aes(x = gDBH, y = n, fill = rank)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(x = gDBH, y = n, label = n),
            position = position_dodge(width = 0.9),
            vjust = -0.5) +
  scale_fill_viridis_d() +
  theme_bw(15) +
  labs(x = "DBH class", y = "Number", fill = "Rank")

Plot BD-2

SpByDBH |> filter(PlotID=="BD-2") |> 
  ggplot(aes(x = gDBH, y = n, fill = rank)) +
  geom_col(position = position_dodge()) +
  geom_text(aes(x = gDBH, y = n, label = n),
            position = position_dodge(width = 0.9),
            vjust = -0.5) +
  scale_fill_viridis_d() +
  theme_bw(15) +
  labs(x = "DBH class", y = "Number", fill = "Rank")

Tree distribution map

BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")

Draw a blank plot with draw_blank_plot_20ha() function.

p <- draw_blank_plot_20ha()
p

Add add tree to map

p + geom_point(data = BDplot.BD1,
               aes(x = gx, y = gy, size = dbh_cm, color =  gDBH),
               alpha = 0.7) +
  xlim(80, 180) + ylim(400, 500) +
  labs(color = "DBH class", size = "DBH(cm)")

For big tree

BDplotA <- BDplot.BD1 |> filter(gDBH == ">=20")
p + geom_point(data = BDplotA, 
               aes(x = gx, y = gy, size = dbh_cm),
               alpha = 0.7, shape = 21,fill = "blue", color = "white") +
    xlim(80, 180) + ylim(400, 500) +
    labs(size = "DBH(cm)")

For medium tree

BDplotB <- BDplot.BD1 |> filter(gDBH == "[5-20)")
p + geom_point(data = BDplotB, 
               aes(x = gx, y = gy, size = dbh_cm),
               alpha = 0.7, shape = 21,fill = "blue", color = "white") +
    xlim(80, 180) + ylim(400, 500) +
    labs(size = "DBH(cm)")

For small tree

BDplotC <- BDplot.BD1 |> filter(gDBH == "[1-5)")
p + geom_point(data = BDplotC, 
               aes(x = gx, y = gy, size = dbh_cm),
               alpha = 0.7, shape = 21,fill = "blue", color = "white") +
    xlim(80, 180) + ylim(400, 500) +
    labs(size = "DBH(cm)")

For a particular species, px: Pinus krempfii

## So do phan bo toan bo cay trong o mau
BDplotSP <- BDplot.BD1 |> filter(scientificName == "Pinus krempfii")

p + geom_point(data = BDplotSP,
               aes(x = gx, y = gy, size = dbh_cm, color =  gDBH),
               alpha = 0.7) +
  xlim(80, 180) + ylim(400, 500) +
  labs(color = "DBH class", size = "DBH(cm)")

We can also try with a detail DBH class

BDplot.BD1 <- BDplot.BD1 |> 
mutate(
  dbh_cm = dbh / 10 , # dbh in mm
  gDBH2 = case_when(
                dbh_cm >= 100 ~ ">=100",
                dbh_cm >= 90 ~ "[90-100)",
                dbh_cm >= 80 ~ "[80-90)",
                dbh_cm >= 70 ~ "[70-80)",
                dbh_cm >= 60 ~ "[60-70)",
                dbh_cm >= 50 ~ "[50-60)",
                dbh_cm >= 40 ~ "[40-50)",
                dbh_cm >= 30 ~ "[30-40)",
                dbh_cm >= 20 ~ "[20-30)",
                dbh_cm >= 10 ~ "[10-20)",
                TRUE ~ "[1-10)")
    )

stemByDBH2 <- BDplot.BD1 |>
  count(gDBH2) |>
  mutate(freq = round(n / sum(n), 3)) |>
  arrange(desc(n))

stemByDBH2
##       gDBH2     n  freq
## 1    [1-10) 11169 0.915
## 2   [10-20)   624 0.051
## 3   [20-30)   246 0.020
## 4   [30-40)   123 0.010
## 5   [40-50)    35 0.003
## 6   [70-80)     3 0.000
## 7     >=100     2 0.000
## 8   [50-60)     2 0.000
## 9   [60-70)     1 0.000
## 10  [80-90)     1 0.000
## 11 [90-100)     1 0.000
stemByDBH2 |>
  ggplot(aes(x = gDBH2, y = n)) +
  geom_col() +
  geom_text(aes(label = n), vjust = -0.5) +
  #scale_fill_manual(values = c("grey","white"))+
  theme_bw(14) +
  labs(x = "DBHClass (cm)", y = "Number of stem")

Dominant species and Biodiversity Index

Dominant species

The Importance Value Index (IVI) in Ecology is the quantitative measure of how dominant a species is in a given ecosystem. IVI is calculated from 3 values:

  1. relative density of each species (RD)

  2. relative basal area of each species (RBA)

  3. relative frequency of each species (this can be omitted in 1 plot)

\[IVI={\left(relative.density + relative.basal.area \right)}\]

Basal area (BA) of each species can be calculated from diameter measured: \(BA = pi * R^2\)

Relative basal area (RBA): \(RBA = BA / sum(BA)\)

Relative density (RD) = number of individuals of each species / total

  BDplot <- BDplot |> 
              mutate(BA = ((dbh/2)/1000)^2*pi) # unit: m2)

BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")
BDplot.BD2 <- BDplot |> filter(PlotID == "BD-2")
IVI.1 <- calcIVI(BDplot.BD1)
head(IVI.1)
## # A tibble: 6 × 7
##   scientificName              n    sG     UT    MD   IVI  rank
##   <chr>                   <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus pierrei       966 6.74  12.7   7.91  20.6      1
## 2 Polyosma nhatrangensis   1175 2.33   4.40  9.63  14.0      2
## 3 Syzygium rubicundum       726 3.21   6.06  5.95  12.0      3
## 4 Pinus krempfii             89 5.28   9.96  0.729 10.7      4
## 5 Castanopsis echinocarpa   539 3.05   5.75  4.42  10.2      5
## 6 Lasianthus saprosmoides  1076 0.236  0.444 8.81   9.26     6
IVI.2 <- calcIVI(BDplot.BD2)
head(IVI.2)
## # A tibble: 6 × 7
##   scientificName              n    sG    UT    MD   IVI  rank
##   <chr>                   <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus coinhensis    357 3.94  7.84   3.08 10.9      1
## 2 Castanopsis echinocarpa   483 3.27  6.51   4.17 10.7      2
## 3 Polyosma nhatrangensis    705 1.31  2.60   6.09  8.69     3
## 4 notdet                    883 0.127 0.253  7.63  7.88     4
## 5 Rhaphiolepis poilanei     436 1.07  2.13   3.77  5.90     5
## 6 Lithocarpus pierrei       262 1.40  2.79   2.26  5.05     6

Figure

figureIVI(IVI.1)

figureIVI(IVI.2)

Diversity Index

biodivIndex(IVI.1)
##   Richness   Simpson  Shannon  Evenness
## 1      124 0.9598452 3.765398 0.7811574
biodivIndex(IVI.2)
##   Richness   Simpson  Shannon  Evenness
## 1      169 0.9756952 4.200136 0.8187561

Biomass and carbon stock

Biomass

  • Above-ground biomass (AGB)

  • Below-ground biomass (BGB)

AGB calculation

Chave et al (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1): 87–99.

  • AGB=Total aboveground biomass (kg) 

  • D=DBH (cm) 

  • H=Height (m) 

  • \(ρ\)=Wood density (g/cm3) 

  • BA=Basal area (cm2) 

Equation
Tropical forest DRY \(AGB = exp(-2.187 + 0.916 * ln(ρ * D^2 * H)\)
\(AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\)
Tropical forest MOIST \(AGB = exp(-2.977 + ln(ρ * D^2 * H))\)
\(AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\) 
library(BIOMASS)  
## For more information on using BIOMASS, visit https://umr-amap.github.io/BIOMASS
## Using temporary cache
##   It is recommended to use a permanent cache to avoid to re-download files on each session.
##   See function createCache() or BIOMASS.cache option.
woodensity <- getWoodDensity(
  genus = BDplot.BD1$genus,
  species = BDplot.BD1$scientificName,
  family = BDplot.BD1$family
)  
## The reference dataset contains 16467 wood density values
## Your taxonomic table contains 124 taxa
head(woodensity)
BDplot.BD1$meanWD <- woodensity$meanWD

\(AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\)

BDplot.BD1 <- BDplot.BD1 %>% 
  mutate(dbh_cm = dbh/10,
         AGB = meanWD * exp(-1.499 + 2.148*log(dbh_cm) + 0.207*log(dbh_cm)^2 - 0.0281*log(dbh_cm)^3)) 
(AGB <- sum(BDplot.BD1$AGB)/1000) #Mg/ha 
## [1] 469.1828

BGB calculation

The ratio developed by Mokany et al. (2006). Critical analysis of root : shoot ratios in terrestrial biomes. Global Change Biology, 12: 84-96 offers specific ratios based on forest type and climate zone and are applicable when the aboveground biomass estimate (shoot) is reported at the stand level

(BGB = 0.275 * AGB)
## [1] 129.0253

Total biomass of 1 hecta

(TB <- AGB + BGB)
## [1] 598.208

Carbon stock

The quantity of carbon can then be estimated by converting biomass to carbon using the IPCC default carbon fraction of 0.47. After that, the carbon stock was multiplied by 44/12 (the carbon atom ratio in the molecular weight of CO2) to get the CO2 equivalent value.

#Carbon stock (ton/ha)
(TC <- TB * 0.47) 
## [1] 281.1578
# CO2 (ton)
(CO2 <- TC * 44/22)
## [1] 562.3155

Carbon credit

1 carbon credit = 1 ton of CO2 equivalent (tCO2e) emissions.

In 2023, Vietnam successfully sold 10.3 million forest carbon credits for the first time through the World Bank at a price of $5/ton (Global average selling price ~$11.2)

1 hecta forest at Bidoup - Nui Ba national park ~ 2800 $

(CO2 * 5)
## [1] 2811.578

Tree distribution map in large scale